1.1 Contributions

Please note that authorship is alphabetical. Contributions are listed below - see github for details and who to blame for what :-).

Ben Anderson (b.anderson@soton.ac.uk @dataknut)

1.2 Code

• Repo: https://github.com/CfSOtago/waterQual/

1.3 Citation

If you wish to refer to any of the material from this report please cite as:

• Anderson, B., (2019) Water Quality in New Zealand: Exploring “official” data , University of Southampton: Southampton, UK.

Report circulation:

• Public.

Report purpose:

• to explore official New Zealand Water Quality data

This work has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 700386 (SPATIALEC).

This work is (c) 2019 the University of Southampton.

2.1 MfE data

Data source: https://data.mfe.govt.nz/data/category/fresh-water/

Data file: mfe-river-water-quality-trends-20082017-19982017-and-19902017-CSV/river-water-quality-trends-20082017-19982017-and-1990-2017-raw-data.csv

df <- paste0(dPath, dataFile)

wqdt <- data.table::fread(df)

Overall there are:

• 1107 sites

This isn’t really trend or raw data - it is a single record per site with indicator of change via a trend slope coefficient or a single % change value. We were looking for the data used to create the time-series plots in (for example) https://www.lawa.org.nz/explore-data/otago-region/river-quality/manuherikia-river/manuherikia-at-ophir/

There are also data quality & missing data issues (see Statistical Annex). It’s noticeable for example that, irrespective of the indicator, over half of stations were thought to have insufficient data to determine trend direction. This was true of 70% of stations for ECOLI and 78% for the Macroinvertebrate Community Index (MCI).

tq <- table(wqdt$TrendCategory, wqdt$npID)

kableExtra::kable(addmargins(tq), caption = "Trend category counts") %>% kable_styling()

kableExtra::kable(100 * prop.table(tq, 2), digits = 2, caption = "Trend category %") %>% 
    kable_styling()

plotDT <- wqdt[, .(meanAnnualSenSlope = mean(AnnualSenSlope, na.rm = TRUE)), 
    keyby = .(npID, TrendCategory)]
ggplot2::ggplot(plotDT, aes(x = TrendCategory, y = meanAnnualSenSlope, group = npID, 
    fill = npID)) + geom_col(position = "dodge")

Figure 2.1: Test data values by indicator

plotDT <- wqdt[, .(meanPercent_annual_change = mean(Percent_annual_change, na.rm = TRUE)), 
    keyby = .(npID, TrendCategory)]
ggplot2::ggplot(plotDT, aes(x = TrendCategory, y = meanPercent_annual_change, 
    group = npID, fill = npID)) + geom_col(position = "dodge")

Figure 2.2: Test data values by indicator

Now try to plot individual trends. This is nasty…

Ideally we want to sort by value within pollutant and ‘fade out’ the ones with insufficient data. Colour coding needs to be careful - for some such as MCI an increase is ‘good’.

p <- ggplot2::ggplot(wqdt, aes(y = sID,
                               x = Percent_annual_change, 
                               alpha = forcats::fct_rev(TrendCategory), # reverse so lightest alpha == least data
                               colour = forcats::fct_rev(TrendCategory)
                               )
                     ) +
  geom_vline(xintercept = 0, colour = "red") +
  geom_point() +
  guides(alpha = guide_legend(title = "Trend and quality code")
         ) +
  guides(colour = guide_legend(title = "Trend and quality code")
         ) +
  scale_color_viridis(discrete = TRUE) +
  theme(legend.position = "bottom") +
  facet_grid(. ~ npID)

p

Figure 2.3: Site level trend plots

4.1 Trend ‘raw’ data

skimr::skim(wqdt)

## Skim summary statistics
##  n obs: 9171 
##  n variables: 42 
## 
## ── Variable type:character ───────────────────────────────────────────────────────────────────────
##        variable missing complete    n min max empty n_unique
##    AnalysisNote       0     9171 9171   2  52     0        9
##         climate       0     9171 9171   2   2     0        6
##   DirectionConf       0     9171 9171   6  22     0       10
##       landcover       0     9171 9171   1   2     0        8
##            npID       0     9171 9171   2   5     0        9
##          Season       0     9171 9171   5   8     0        3
##             sID       0     9171 9171   8  11     0     1107
##      src_of_flw       0     9171 9171   1   2     0        5
##   TrendCategory       0     9171 9171  10  17     0        3
##  TrendDirection       0     9171 9171  10  13     0        4
## 
## ── Variable type:integer ─────────────────────────────────────────────────────────────────────────
##         variable missing complete    n       mean         sd     p0
##                D     123     9048 9171    3834.87    8388.81     15
##  no_censorlevels       0     9171 9171       0.64       1.05      0
##             nObs       0     9171 9171     119.47      89.47      5
##          nzreach       0     9171 9171 8353123.65 4569119.66  1e+06
##        nzsegment       0     9171 9171 8463644.58 4596199.64  1e+06
##           Period       0     9171 9171      13.97       6.03     10
##                S     123     9048 9171    -176.53    1953.93 -24882
##                X       0     9171 9171    1947.94    1673.87      1
##      p25     p50        p75        p100     hist
##      463     703 4186       56280       ▇▁▁▁▁▁▁▁
##        0       0    1          21       ▇▁▁▁▁▁▁▁
##       40     117  120         336       ▆▁▇▁▁▃▁▁
##    3e+06   8e+06    1.3e+07     1.5e+07 ▃▇▂▆▂▂▆▆
##  3149001 8206040    1.3e+07     1.5e+07 ▃▇▂▃▅▂▇▆
##       10      10   20          28       ▇▁▁▁▃▁▁▁
##     -133     -14   50       26493       ▁▁▁▇▁▁▁▁
##      473    1447 3290.5      5583       ▇▃▃▂▂▂▂▂
## 
## ── Variable type:logical ─────────────────────────────────────────────────────────────────────────
##  variable missing complete    n mean                       count
##  Seasonal       0     9171 9171 0.58 TRU: 5278, FAL: 3893, NA: 0
## 
## ── Variable type:numeric ─────────────────────────────────────────────────────────────────────────
##               variable missing complete    n        mean        sd
##         AnnualSenSlope     124     9047 9171      -0.16       5.78
##              Intercept     124     9047 9171      36.97     121.26
##                    lat       0     9171 9171     -40.86       3.09
##                    Lci     124     9047 9171      -1.58       9.49
##                   long       0     9171 9171     173.75       2.64
##                 Median     123     9048 9171      35.61     113.73
##                  nzmge       0     9171 9171 2581624.68  218056.88
##                  nzmgn       0     9171 9171   6e+06     342744.62
##                  nztme       0     9171 9171 1672489.23  219604.57
##                  nztmn       0     9171 9171 5473697.5   344920.76
##                      p     123     9048 9171       0.29       0.32
##  Percent_annual_change     124     9047 9171      -0.65       4.25
##            Probability     123     9048 9171       0.57       0.38
##         Probabilitymax     124     9047 9171       0.62       0.38
##         Probabilitymin     124     9047 9171       0.5        0.39
##          prop_censored       0     9171 9171       0.076      0.18
##            prop_unique       0     9171 9171       0.58       0.25
##        Sen_Probability     124     9047 9171       0.57       0.39
##               Sen_VarS     124     9047 9171  147300.04  533409.77
##                    tau     123     9048 9171      -0.031      0.19
##                    Uci     124     9047 9171       1.1        7.77
##                   VarS     123     9048 9171  146978.2   533006.67
##                      Z     123     9048 9171      -0.38       2.75
##             p0          p25           p50           p75       p100
##     -201.49         -0.005       -5.9e-06       0.0023      135.16
##     -333.91          0.02         0.36          3.89       2416.87
##      -46.63        -43.61       -40.26        -38.29        -35.04
##     -317.46         -0.076       -0.0017       -2.5e-05      62.14
##      167.45        172.04       174.66        175.69        178.43
##    5e-04             0.02         0.38          4.04       3000   
##  2080600       2427919      2658279       2746800         3e+06   
##  5390200       5733145.39   6100300       6316100       6684603   
##  1175003       1519856.81   1749951.5     1836751.23    2078090   
##  4828129.84    5170669.03   5538610.78    5761606.95    6122361.21
##    3e-68             0.0084       0.16          0.52          1   
##      -80.99         -1.82        -0.025         0.87         24.91
##        1.5e-68       0.16         0.65          0.96          1   
##    1e-67             0.23         0.78          0.99          1   
##        0             0.063        0.53          0.92          1   
##        0             0            0             0.03          0.98
##        0.0083        0.4          0.62          0.78          1   
##        1.5e-68       0.16         0.67          0.97          1   
##       28.33       1267.17      4221.67      30193.83    4233425.33
##       -0.79         -0.15        -0.033         0.08          0.8 
##      -99.32          0            0.001         0.062       242.93
##       28.33       1245.33      4183.17      30167.5     4233425.33
##      -13.39         -1.77        -0.39          1            17.46
##      hist
##  ▁▁▁▁▇▁▁▁
##  ▇▂▁▁▁▁▁▁
##  ▅▅▅▅▆▇▆▂
##  ▁▁▁▁▁▁▇▁
##  ▂▂▃▃▃▇▆▁
##  ▇▁▁▁▁▁▁▁
##  ▂▂▂▃▃▇▃▁
##  ▅▃▅▅▆▇▅▂
##  ▂▂▂▃▃▇▃▁
##  ▅▅▅▅▇▇▆▂
##  ▇▂▂▁▁▁▁▁
##  ▁▁▁▁▁▃▇▁
##  ▅▂▁▂▂▂▂▇
##  ▃▁▁▁▁▁▂▇
##  ▇▂▂▂▂▂▂▇
##  ▇▁▁▁▁▁▁▁
##  ▂▃▃▅▆▇▇▅
##  ▅▂▁▁▁▂▂▇
##  ▇▁▁▁▁▁▁▁
##  ▁▁▃▇▆▂▁▁
##  ▁▁▇▁▁▁▁▁
##  ▇▁▁▁▁▁▁▁
##  ▁▁▃▇▂▁▁▁

tn <- wqdt[!is.na(AnalysisNote) & AnalysisNote != "ok", .(sID, npID, AnalysisNote)]

kableExtra::kable(tn, caption = "Analysis notes")